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ABSTRACT 

We study the photometric properties of Lyman-break galaxies (LBGs) formed by 
redshift z — 3 in a set of large cosmological smoothed-particle hydrodynamics (SPH) 
simulations of the A cold dark matter (ACDM) model. Our numerical simulations 
include radiative cooling and heating with a uniform UV background, star formation, 
supernova feedback, and a phenomenological model for galactic winds. Analysing a 
series of simulations of varying box size and particle number allows us to isolate the 
impact of numerical resolution on our results. We compute spectra of simulated galax- 
ies using a population synthesis model, and derive colours and luminosity functions 
of galaxies at z — 3 after applying local dust extinction and absorption by the in- 
tergalactic medium (IGM). We find that the simulated galaxies have U n — G and 
G — R colours consistent with observations, provided that intervening absorption by 
the IGM is applied. The observed properties of LBGs, including their number density, 
colours, and luminosity functions, can be explained if LBGs are identified with the 
most massive galaxies at z = 3, having typical stellar mass of M* ~ 10 10 /i _1 Mq, a 
conclusion broadly consistent with earlier studies based on hydrodynamic simulations 
of the ACDM model. We also find that most simulated LBGs were continuously form- 
ing stars at a high rate for more than one Gyr up until z = 3, but with numerous 
starbursts lying on top of the continuous component. Interestingly, our simulations 
suggest that more than 50% of the total stellar mass and star formation rate in the 
Universe are accounted for by galaxies that are not detected in the current generation 
of LBG surveys. 

Key words: cosmology: theory - galaxies: evolution - galaxies: formation - methods: 
numerical. 



1 INTRODUCTION 

Based on surveys made with modern large telescopes, 
the number of observed Lyman-br eak galaxies (LBGs; e.g. 
ISteideL Petting fc Ham ilton! Il995l) has reached about ~ 
1000 (Ste idel et al J 12003') by now, providing a unique win- 
dow for studying early galaxy formation. In these obser- 
vations, a large sample of LBG-candidates is first selected 
as "drop-outs" in certain sets of colour-colour planes (e.g. 
Un — G vs. G — R), and then later spectroscopic follow-up 
observations are performed to obtain the redshifts of the 
identified candidates. Most surveys to date have been car- 
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ried out at optical wavelengths, corresponding to the rest- 
frame far ultra-violet (UV) at redshift z = 3. The L BGs 
found in this way show strong clustering (lAdelberger et alJ 
119981: iGiavalisco et alJll99S! : ISteidel et al]|l99S!h which has 
generally been interpreted as indirect evidence that LBGs 
reside in massive dark matter halos. Several studies of 
LBGs based on semi-analytic models for galaxy formation 
agree with the massive dark matter halo hypothesis (e.g. 
| Mo fc Fukugitalll996t feaugh et alJll998t iKauffmann et alJ 
Il999t iMo. Mao, fc Whitelll999lh but a numerical study by 
IJing fc SutoT " lH99c1) ~using collisionless N-body simulations 
expressed some concerns whether the clustering of dark mat- 
ter halos of < 10 12 /i -1 Mq can explain the observed clus- 
tering strength of LBGs. However, a later investigation by 
iKatz. Hernauist fc Weinberg! Jl999f) argued that the CDM 
models have no difficulty in explaining the strong observed 
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clustering of LBGs, although the substantial box size effects 
in their analysis introduced some uncertainty. 

More recently, LBGs have also been stud- 
ied in the near infrared (IR~) llSawicki fc Yed 199S ; 



Parjovich^Dickinson, & Ferguson 200 1|: [ Pettini et aO 



2001; 



Shaplev et alJl200lt iRudnick et alJl200ll) . Near IR obser- 
vations are less affected by dust obscuration and directly 
probe the rest-frame optical properties of LBGs at z = 3. 
Therefore, they allow a more reliable derivation of the 
stellar masses of LBGs. Based on these observations, it has 
been suggested that the value of the extinction lies in the 
range E(B — V) » 0.0 — 0.3 with a median value of 0.15, 
and that a significant old stellar population exists in L BGs 
at z = 3, with a stellar m ass of w 1O 1O M iPapovich et alj 
1200 ll: IShaplev et all 1200 j) . The existence of such a stellar 
component would again suggest that LBGs are embedded in 
massive dark matter halos which have continuously formed 
stars over an extended period of roughly one Gyr up to 
z = 3. These systems would then most likely evolve into 
elliptical galaxies at the present day, or into the spheroidal 
components of massive spiral galaxies. 

However, in a competing model, LBGs have been sug- 
gested to be merger-i nduced starbursting systems associated 
with l ow-mass halos (lL^wentJial^t_a lJ[l997|; ISawicki fc Yed 
Il998l: ISomerville. Primack. fc Faberl 2001), and in some 
cases merger-induced starbursts are given a crucial role even 
when LBGs are the m ost massive galaxies at their time 
ISomerville et alJl200lfl . Given that the merger rate is ex- 
pected to be quite high at z ~ 3, these scenarios provide 
an interesting alternative to the more conventional picture 
which associates LBGs with the most massive systems. 

Self-consistent hydrodynamic simulations are an ideal 
tool for trying to distinguish bet ween these d i fferen t sce- 
narios for the nature of LBGs. [ Dave et alJ lll999l) and 
IWeinberg. Hernquist fc Katd J2002T) were the first to em- 
ploy smoothed particle hydrodynamics (SPH) simulations 
to this end. However, their box size of 11.1 ft -1 Mpc repre- 
sented an important limitation, because the space-density 
of LBGs is so low that only a few of them can be found in 
a volume of this size, as we will discuss further in Section 0] 
Therefore, simulations with a larger box si ze are desirable to 
obtain larger samples of simulated LBGs. iNagamind' l|2002t) 
used an Eulerian hydrodynamic simulation with a box size 
of Lbox = 25 ft -1 Mpc, tracing the merger history of galaxies 
from z = 5 to z = 0. The results of these earlier numerical 
studies were consistent with each other, and agreed reason- 
ably well with the observations, within the uncertainties. 
In particular, the median stellar masses of LBGs were pre- 
dicted to be ~ 10 10 /i _1 Mq, and the simulated galaxies were 
experiencing significant star formation rates (> 30MQ/yr) 
for extended periods of time (> 1 Gyr). 

In this paper, we improve on the earlier numerical stud- 
ies of LBGs by using a new set of high-resolution numerical 
simulations. These simulations are based on a novel model 
for the physics of star formation and feedback, and they 
use a more accurate implementation of SPH. For the first 
time, we also systematically study the effects of resolution 
and box size in the context of simulated LBG galaxies. The 
treatment of star formation and feedback we use is based on 
a sub-resolution multi-phase description of the dense, star- 
forming interstellar medium (ISM), and a phenomenologi- 
cal model for strong feedback by galactic winds, as recently 



proposed bv lSpringel fc Hernauistl (l2003aTl . This model has 
been shown to provide converged star formation rates for 
well-resolved galaxies, with a cosmic star formation history 
consistent with recent observatio ns llSpringel fc Hernauistl 
l2003btlHernauist fc Springelll2003t) . 3 The inclusion of winds 
was motivated by the realisation that galactic outflows at 
high redshift dPettini et all [20021 likely play a key role 
in distributing meta ls into the intergalactic medium (e.g. 
lAguirre et al]l200ldlbh . as well as being important for the 
regulation of star formation activity. In fact, winds may 
also alter the distribut ion of neutral gas around galaxies 
l|Adelb crger et al.ll2003|) . although the d etails of how this 
process may happen r emain u nclear [e.g. ICroft et al.ll2002l : 
iKollmeier et alJ 2003; Br^gcoJj^taL|]2003J). I n passing, we 
note t hat both iDesiacaues et alJ l|2003h and iMaselli et alJ 
( 2003) have found that t he Lyman-q transmissiv ity close 
to LBGs, as measured bv lAdelberger et alJ (I2003F) . is bet- 
ter reproduced if LBGs are i dentified as dwa r f starb urst- 
ing galaxies as proposed in ISomerville et alJ (1200 if) and 
Weathcrlev & Warren (2003). We will discuss the work of 



Weatherlev fc WarrerJ (|2CJ03) in Section 03 Together with 



the increase in numerical resolution provided by our simu- 
lations, it is of interest to see how our refined physical mod- 
elling modifies the predictions for LBG properties within the 
ACDM scenario. 

This paper is organised as follows. In Section [5] we 
briefly introduce the numerical parameters of our simula- 
tion set. In Section^ we then describe our method for com- 
puting spectra of simulated galaxies both in the rest-frame 
and the observed frame. In Sections |1] and 03 we show the 
colour-colour diagrams and colour-magnitude diagrams of 
simulated galaxies, and we discuss the number density of 
colour-selected LBGs, as well as the stellar masses of LBGs 
at z = 3. We then investigate the rest-frame l/-band lu- 
minosity function and observed i?-band luminosity function 
in Section HJ followed by an analysis of the star formation 
histories of LBGs in Section |7| Finally, we summarise and 
discuss the implications of our work in Section [8] 



2 SIMULATIONS 

We analyse a large set of cosmological SPH simulations with 
varying box size, mass resolution and feedback strength, 
as summarised in Table Our box size ranges from 10 
to 100 ft -1 Mpc on a side, with particle numbers between 
2 x 144 3 and 2 x 486 3 , giving gaseous mass resolutions 
in the range 3.3 x 10 5 to 3.3 x 10 8 /i _1 M Q . These simu- 
lations are partly tak en from a study of the cosm ic star 
formation history by ISpringel fc Hernauistl ll200 3b'l. sup- 
plemented by additional runs with weaker or no galac- 
tic winds. A similar set of simulatio ns was used by 
iNagamine. Springel. fc Hernauistl d2004 a b) to study the 
properties of damped Lyman-a absorbers, but here we anal- 
yse the 'G6'-run which has higher resolution than the 'G4'- 



3 We note an error in figure 12 of lSpringel fc Hernauistl <2003tft 
in which the observational estimates of the SFR were plotted too 
high by a factor of h" 1 = 1.4. When corrected, the observed 
points are in better agreement with the theoretical estimates; see 
astro-ph/0206395 
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Run 


Box Size 


V 

"p 


m DM 


' ' 'saw 




■^end 




03 


10.00 


2 x 144 3 


2.42 X 10 7 


3.72 x 10 6 


2.78 


2.75 


none 


P3 


10.00 


2 X 144 3 


2.42 x 10 7 


3.72 x 10 6 


2.78 


2.75 


weak 


Q3 


10.00 


2 X 144 3 


2.42 X 10 7 


3.72 x 10 6 


2.78 


2.75 


strong 


Q4 


10.00 


2 X 216 3 


7.16 x 10 6 


1.10 x 10 6 


1.85 


2.75 


strong 


Q5 


10.00 


2 X 324 3 


2.12 x 10 6 


3.26 x 10 5 


1.23 


2.75 


strong 


D4 


33.75 


2 X 216 3 


2.75 x 10 8 


4.24 X 10 7 


6.25 


1.00 


strong 


D5 


33.75 


2 X 324 3 


8.15 x 10 7 


1.26 x 10 7 


4.17 


1.00 


strong 


G5 


100.0 


2 x 324 3 


2.12 x 10 9 


3.26 x 10 s 


8.00 


0.00 


strong 


G6 


100.0 


2 x 486 3 


6.29 x 10 s 


9.67 x 10 7 


5.00 


0.00 


strong 



Table 1. Simulations employed in this study. The box size is given in units of h~ 1 Mpc, N p is the particle number of dark matter and 
gas (hence X 2), muM and "igas are the masses of dark matter and gas particles in units of /i _1 M0, respectively, e is the comoving 
gravitational softening length in units of h~ 1 kpc, and z on( j is the ending redshift of the simulation. The value of e is a measure of spatial 
resolution. From the top to the bottom row, we collectively call the first 5 runs (03 to Q5) 'Q-series', D4 & D5 are called 'D-series', and 
G5 & G6 are called 'G-series'. 



run used in the previous studies. The simulations with the 
same box size are run with the same initial condition. 

There are three main novel features to our simulations. 
First, we use the new "conservati ve entropy" formulation of 
SPH JSprineel fc Hernauistj2 002) which explicitly conserves 
entropy (in regions without shocks), as well as momentum 
and energy, even when one allows for ful ly adaptive smooth- 
ing lengths (see e.g. iHernauistl Il993l) . This formulation 
moderates the overcoolin g problem present in earlier formu- 
lations of SPH (see al so lYoshida et alJl2002l : iPearce et alJ 
ll999l:ICroft et alj|200ll) . 

Second, highly over-dense gas particles are treated with 
an effecti ve sub-resolution multiphase ISM model, as de- 
scribed bv lSpringel fc Herncmistl J2003aft . In this model, each 
gas particle represents a statistical mixture of cold clouds 
and a hot ambient phase. Cold clouds can grow by radia- 
tive cooling out of the hot medium, and form the reservoir 
of neutral gas for star formation. Once star formation takes 
place, supernova explosions deposit energy into the hot gas, 
heating the gas and evaporating the cold clouds, transferring 
cold gas back into the hot phase. This feedback establishes 
a self-regulated cycle of star formation. 

Third, a phenomenological model for galactic winds is 
implemented. In this model, gas particles are stochastically 
driven out of the dense star-forming region by assigning an 
extra momentum in random directions, with a rate and mag- 
nitude chosen to reproduc e mass-loads and wind speeds s im- 
ilar to those observed. (See lSpringel fc Hernauistl (l2003aT) for 
a detailed discussion of this method.) Most of our simula- 
tions employ a "strong" wind of speed 484 km s -1 , but for 
the 03, P3, & Q3-runs, we also varied the wind strength to 
examine the effect of feedback from galactic winds. The runs 
with 10 h~ Mpc boxes are collectively called the 'Q-series', 
and the resolution is increased from Q3 to Q4, and to Q5. 
The other series ('D-'and 'G-Series') extend the strong wind 
results to larger box sizes and hence lower redshift. Our nam- 
ing convention is such that runs of the same model (box size 
and included physics) are designated with the same letter, 
with an additional number specifying the particle resolution. 

Our calculations include a uniform UV back g round 
radiation field with a modified lHaardt fc Madaul (Il996l) 
spectrum, where reionisation takes place at z ~ 6 (see 
iDave et al]|l999Fl as suggested by the quasar observations 



(e.g. iBecker et al"ll200ll) and radiative transfer calculations 
of the impact of the stellar source s in our simulations on 
the IGM (e.g. ISokasian et al.lEoO^) . The early reionisation 
at higher redshift, as suggested by the WMAP satellite, 
should not affect the results presented in this paper because 
we are mainly dealing with halos with virial temperatures 
above 10 4 K and the infalling gas will radiatively cool even 
if it was photoionised by th e reionisation at z > 10 be - 
fore falling into the halo (see ISpringel fc Hernauistll2003b() . 
The radiat ive cooling and he ating rate is computed as de- 
scribed by iKatz et alJ il996T) assuming that the gas is op- 
tically thin and in ionisation equilibrium. The adopted cos- 
mological parameters of all runs are (f2 m , Qa, fit>, ""8, h) = 
(0.3,0.7,0.04,0.9,0.7). The simulations were performed on 
the Athlon-MP cluster at the Center for Parallel Astrophys- 
ical Computing (CPAC) at the Harvard-Smithsonian Cen- 
ter for Astroph ysics, using a modified version of t he parallel 
GADGET code iSpringel. Yoshida fc Whitdl200ll) . 



3 METHOD 

We identify simulated galaxies as isolated groups of stars 
using a simpli fied varia n t of t he SUBFIND algorithm pro- 
posed bv lSpringel et alJ (2001). In detail, we first compute 
an adaptively smoothed baryonic density field for all star 
and gas particles, allowing us to robustly identify centres 
of individual galaxies as isolated density peaks. We find the 
full extent of these galaxies by processing the gas and star 
particles in the order of declining density, adding particles 
one by one to the galaxies. If all of the 32 nearest neighbours 
of a particle have lower density, this particle is considered to 
be a new galaxy 'seed'. Otherwise, the particle is attached to 
the galaxy that its nearest denser neighbour already belongs 
to. If the two nearest denser neighbours belong to different 
galaxies, and one of these galaxies has less than 32 particles, 
these galaxies are merged. If the two nearest denser neigh- 
bours belong to different galaxies and both of these galaxies 
have more than 32 particles, then the particle is attached 
to the larger group of the two, leaving the other one in- 
tact. Finally, we restrict the set of gas particles processed 
in this way to those particles which have at least a density 
of 0.01 times the threshold density for star formation, i.e. 
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Pth = 8.6 x 10 6 /i 2 M Q kpc- 3 fsee lSpringel fc Hernauistl2003al. 
for details on how this parameter is determined). Note that 
we are not interested in the gaseous components of galax- 
ies - we only include gas particles because they make the 
method more robust. Since most galaxies contain very dense 
star-forming gas, such a method make it particularly easy 
to select galaxies when the gas is included. 

We found that the above method robustly links up star 
particles that belong to the same isolated galaxy. A simpler 
FoF algorithm with a small linking length can achieve a 
largely similar result, but the particular choice for the link- 
ing length one needs to make in this method represents a 
problematic compromise, either leading to artificial merging 
of galaxies if it is selected too large, or to loss of star parti- 
cles that went astray from the dense galactic core, if selected 
too small. Note that, unlike in the detection of dark mat- 
ter substructures, no gravitational unbinding algorithm is 
needed to define the groups of stars that make up the galax- 
ies formed in the simulations. 

We only consider galaxies with at least 32 particles (star 
and gas combined) in our subsequent analysis. Each stellar 
particle contained in them is tagged by the simulation code 
with its mass, formation time, and metallicity of the gas 
particle that it formed out of. Based on these three tags, we 
compute the emission from each stellar particle, and co-add 
the flux from all particles for a given galaxy to obtain the 
spectrum of the simulated galaxy. We use the population 
synthesis mo del GISSEL99 llBruzual fc Chariolll993tl that 
assumes the ISalpete 3 Jl995l) initial mass function with a 
mass range of [0.1, 125] Mq. 

Once the intrinsic s pectrum is compute d, we apply the 
Calzetti extinction law | Calzettietal . 2000) with three dif- 
ferent values of E(B — V) = 0.0,0.15,0.3 to investigate 
the effects of extinction. These values span the range of 
E( B — V) estimated from observations of LBGs at z = 
3 iAdelbereer fc Steidell bood IShaplev et al.ll200ll) . Rest- 
frame colours and luminosity functions of the simulated 
galaxies are then derived using the spectra computed in this 
manner. 

To obtain the spectra in the observed frame, we red- 
shift the spectra and appl y absorption b y the IGM fol- 
lowing the prescription by iMadaul il995T) . Once the red- 
shifted spectra in the observed frame are obtained, we con- 
volve them with various fil ter functions, including U n ,G,R 
JSteidel fc Hamilton! Il993h and standard Johnson bands, 
and compute the magnitudes in both AB and Vega systems. 
Apparent U n , G, R magnitudes a re computed in the A B sys- 
tem to compare our results with lSteidel et all j2003t) . while 
the rest-frame V-band magnitude is comput ed in the Vega 
system to compare our results with those of IShaplev et alJ 
(2001). 



4 COLOUR-COLOUR DIAGRAMS 

In Figure we show the colour-colour magnitude diagrams 
of simulated galaxies at z — 3 in the U„ — G vs. G — R 
plane, for ga laxies brighter than R = 25.5 (the magnitude 
cut used by ISteidel et al.l I2003T) . The three different sym- 
bols represent three different values of Calzetti extinction: 
E(B - V) = 0.0 (blue dots), 0.15 (green crosses), and 0.3 
(red open squares). The long-dashed lines mark the colour- 



colour selection criteria applied by ISteidel et alJ ll2003f) to 
identify LBG candidates at z ~ 3. 

It is encouraging to see that in all the panels most of the 
simulated galaxies actually satisfy the observational colour 
selection criteria. This suggests that the simulated galaxies 
have realistic colours compared to the observed LBGs at 
z = 3. In runs with larger box size (D- and G-series), the 
distribution is wider than the one in the Lbox = 10 ft" 1 Mpc 
runs (Q-series), and the distribution actually extends be - 
yond the colour selection boundaries of lSteidel et alJ {2003). 
The box size of the Q-series (Lbox = 10 /i -1 Mpc) is too 
small to contain a significant sample of LBGs, and no LBGs 
with R < 25.5 can be found in the Q-series for the cases of 
E(B - V) = 0.15 and 0.3. 

If we had not applied the Madau (1995) absorption to 
the spectra, the distribution would not have fallen into the 
colour selection region. This is because the simulated galax- 
ies would then appear too bright in the UV, such that the 
distribution would fall below U n — G < 1.0. As the level of 
extinction by dust is increased from E(B — V)= 0.0 to 0.3, 
the measured points move towards the upper right corner of 
each panel. This behaviour is expected for a standard star- 
forming galaxy spe ctrum, as demonstrated in Figure 2 of 
ISteidel et all J20031) . At the same time, the number of galax- 
ies that satisfy the magnitude cut of R = 25.5 decreases with 
increasing extinction, because stronger extinction results in 
a redder galaxy spectrum. 

Note that increasing the resolution from Q3 to Q4, and 
then to Q5, also reduces the number of galaxies brighter than 
R — 25.5 slightly. The same is true for D4 and D5 runs. As 
the resolution is increased, not only are the most massive 
halos in the simulation better resolved, but also all of their 
progenitors. The better resolution then allows a more ac- 
curate treatment of the wind-feedback in these progenitor 
generations of galaxies; the net result of this is a decrease in 
the final luminosity of the brightest galaxies. For the 'Q5'- 
run, there are actually no galaxies brighter than R — 25.5. 

This is not the case for the lower resolution, larger box 
size simulations where the trend appears to reverse. In the 
case of the G-series, the number of galaxies brighter than 
R = 25.5 actually increases as the resolution increases. As 
we will show in Section [6.21 this is because the peak of the 
simulated i?-band luminosity function is still on the brighter 
side of R = 25.5, and the increase in the number of galaxies 
near the peak of the luminosity function wins over the slight 
decrease of the number at the brightest end. 

We now investigate the number density of LBGs in 
the simulation. In Figure |5J we plot the number density 
of galaxies that satisf y the colour-colour selection criteria 
of lSteidel etai] J2003D . for all the runs shown in Figure 
Three different symbols represent the three different val- 
ues of extinction we used: E(B — V) = 0.0 (black open 
squares), 0.15 (blue filled squares), and 0.3 (red filled tri- 
angles). The points for the same value of E(B — V) are 
connected to guide the eye. All the cases with zero LBGs 
are indicated as N = —3.8 dex for plotting purposes. A 
conservative range for the observed number densities of 
LBGs is shown as a shaded region, with a median value of 
4 x 10~ 3 fo 3 Mpc~ 3 fAdelberger et al.ll2003t) . We note that 
iGiavalisco fc Dickinson! i200ll) reported a slightly smaller 
value of 2 - 3 x 10" 3 /i 3 Mpc" 3 . 

From this figure, we find that the D4 & D5 runs with 
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Figure 1. Colour-colour diagrams (U n — G vs. G — R) of simulated galaxies at z = 3 with apparent R magnitudes brighter than R < 25.5. 
Three different symbols represent three different values of extinction: E(B — V) = 0.0 (blue dots), 0.f5 (green crosses), and 0.3 (red open 
squares). The number of simulated galaxies of e ach type is indicate d from top to bottom in each panel, respectively. The dashed lines 
mark the colour-colour selection criteria used bv lSteidel et alj J2003fl to identify LBG candidates at z ~ 3. 



E{B - V) = and the G5 & G6 run with E(B -V) = 0.15 
appear to have reasonable number densities, while all runs 
with E(B — V) — 0.3 seem to underpredict the number 
density. The 03 (no-wind) and P3 (weak wind) runs with 
E(B — V) = overpredict the number density significantly. 
Apparently, without strong feedback by galactic winds and 
extinction, the LBGs simply become too bright and hence 
too abundant above a given brightness limit. Other null re- 
sults of the Q-series are affected by the small box size of the 
simulation. We will show the effect of the box size more ex- 
plicitly when we discuss the luminosity function of galaxies 
in Section [S] 



5 COLOUR-MAGNITUDE RELATION AND 
STELLAR MASSES 

In Figure we show the distribution of galaxies at z — 3 
in the 'G6'-run on the plane of J?-band apparent magnitude 
and G—R colour. We here chose the 'G6'-run because it gives 
a reasonable fit to the observed luminosity functions and the 
number density, and it has higher resolution than the G5 



run. Again, we use three different symbols for three different 
values of extinction: E(B — V) = 0.0 (blue dots), 0.15 (green 
crosses), and 0.3 (red open squares). The long-dashed lines 
an d the arrows indicat e the colour-selection criteria applied 
bv lSteidel et alJ i2003tl to select LBG candidates at z ~ 3. 

We see that most of the galaxies brighter than R = 25.5 
automatically satisfy the criterion G — R < 1.2, and only a 
small fraction of galaxies with R < 25.5 fall out of the region. 
There is a significant population of dim (7? > 28) galaxies 
with G — R > 1.2. We will see below that these are low-mass 
galaxies with stellar masses M* < 10 10 /i _1 M Q . 

In Figure 0] we show the J?-band apparent magnitude 
vs. stellar mass of simulated galaxies at z = 3. As before, 
we plot results for three values of extinction, using different 
sy mbols. We also mar k the magnitude limit R < 25.5 used 
bv lSteidel et alJ <l2003l) with a vertical long-dashed line and 
arrows. 

From this Figure, we see that the LBGs in the Q- 
series with R < 25.5 have typically stellar masses in the 
range M* = 10 9 - 10 10 ft _1 M , while those in 'D5' take a 
somewhat wider interval, nearly covering the range 10 8 — 
10 11 /i _1 Mq. The G6 run covers the even wider mass range 
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Figure 2. Number density of galaxies for all the runs shown in 
Figure Q Three different symbols denote three different values of 
extinction: E(B — V) = 0.0 (black open square), 0.15 (blue filled 
square), and 0.3 (red filled triangle). Points with the same value of 
E(B — V) are connected to guide the eye. A conservative range of 
the observed number density of LBGs is shown as the shaded re- 
gion, with a central value of 4 X 10 — 3 /i 3 Mpc — 3 iAdclbcrge r et all 
2003). All the cases with zero LBGs are indicated as N = —3.8 
dex for plotting purposes. 
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Figure 3. Colour-magnitude diagram for the 'G6'-run at z = 3. 
The distribution of galaxies is shown in the plane of R apparent 
magnitude and G — R colour. Three different symbols give results 
for three different values of extinction: E(B — V) = 0.0 (blue 
dots), 0.15 (green crosses), and 0.3 (red open squares). The long- 
dashed lines and the arrow indicate the colour-selection criteria 
used bv lSteidel et alJ <2003ft to identify LBG candidates at z ~ 3. 



-Mstar = 10 s — 10 12 /i _1 Mq depending on the value of extinc- 
tion. As the wind strength is increased from 03 to P3, and 
then to Q3-run, we observe that the distribution of galaxies 
becomes slightly sparser at the most luminous (i.e. massive) 
end of the distribution. We see the same effect more clearly 
in the luminosity functions discussed in Section HJ Note that 
since the numerical resolution is identical for 03, P3, and 
Q3, the total baryonic mass in the simulation box is the same 
for these runs, but the runs with stronger winds convert less 
of their baryonic mass components into stars. 

In simulations with larger box size (D- and G-series), 
the masses of the most luminous galaxies found are sub- 
stantially larger compared with those in the Q-series. This 
is primarily a result of a finite box size effect; the brightest 
galaxies have a very low space-density, such that they are 
simply not found in simulations with too small a volume. 
On the other hand, the mass resolution of the small box 
size simulations is much better, so that they can resolve a 
much larger number of dim low-mass galaxies. Note that the 
lower cut-off of the distribution in each panel is given by the 
stellar particle mass of the corresponding simulation. 

Interestingly, we find that the amount of stellar mass 
contained in the LBGs that satisfy the colour selection cri- 
teria is in general far from being close to the total stellar 
mass. In the 'D5'-run, for example, the LBGs that pass the 
criteria contain only 20% of all the stars for an extinction of 
E(B-V) = 0.0. In the 'G5'-run, the fraction is 57% and 20% 
for E(B - V) = 0.0 and 0.15, respectively. In the 'G6'-run, 
the fraction is 52% and 27% for E(B - V) = 0.0 and 0.15, 
respectively. This means that the current surveys of LBGs 
fail to ac count for more than h alf of the stellar mass in the 
Universe. iNagamine et alJ (120041) reach a similar conclusion 
from different theoretical arguments by comparing the re- 
sults of two different type s of hydrodynamic simulation s and 
the theoretical model of iHern quist fc Springell J2003I) with 
near-IR observation s of galaxies (e.g. lDickmsonet^n1l2CJ0"3l 
iRudnick et"ai]l2003l) . Such missed stellar masses might be 
hidden in the red population of galaxies as suggested by 
iFranx et all §003), but in our simulations the missed stel- 
lar mass fractions are mainly contained in the faint galaxies 
below the observational magnitude limit. 



6 LUMINOSITY FUNCTIONS 

6.1 Rest- frame V-band luminosity function 

In Figure|SJ we show the rest-frame V-band luminosity func- 
tions for all the simulated galaxies at z = 3. We have also 
tried restricting the luminosity functions to only those galax- 
ies that satisfy the colour-selection criteria. Most of the 
galaxies that fall out of the colour-selection criteria are the 
less luminous ones as demonstrated in Figure This only 
affects the faint-end of the luminosity function not probed 
by current LBG surveys. Therefore we show the luminosity 
functions for all galaxies in our simulations in Figures Qj] |S| 
&Qfor clarity. As before, we plot results for three different 
values of extinction: E(B — V) = 0.0 (blue long-dashed), 
0.15 (green solid), and 0.3 (red dot-dashed). Note that 
stronger extinction simply shifts the luminosity function to 
the fainter side, without changing its shape. The black short- 
dashed line included in the figure gives the observational 
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Figure 4. Stellar mass vs. apparent magnitudes of simulated galaxies. Three different symbols are used to show results for three different 
values of extinction: E(B — V) = 0.0 (blue dots), 0.15 (g r een cr osses), and 0.3 (red open squares). The vertical long-dashed line and the 
arrow indicates the magnitude cut used bv lSteidel 



estimate of the rest -frame V-band luminosity function by 
IShanlev et alJ (j200lfl with a faint-end slope of a = —1.85, 
normalisation $* = 0.18 x 10~ 2 /i 3 Mpc -3 , and characteristic 
magnitude M* = -22.21 + 5 log h = -22.98, for h = 0.7. 
The observational magnitude limit of Mv = —20.5 is shown 
by the arrows. 

The most prominent feature seen in all panels is that the 
luminosity functions of the simulated galaxies are all very 
steep, with a faint-end slope comparable to a ~ —2, which 
is the slope of the dark matter halo mass function. This sug- 
gests that the strong feedback included in the simulations 
has not been able to reduce the luminosities of low-mass 
galaxies much more strongly than those of more massive 
systems; if such a differential effect existed, it should have 
manifested itself as a flattening of the faint-end compared to 
the halo mass function. At face value, however, the observa- 
tional data actually support a rather steep faint-end slope 
at z — 3, quite close to that of the halo mass function. One 
should note, however, that the observational estimate of the 
slope a at z — 3 is very uncertain because the observations 
can only reach down to a magnitude of Mv ~ —20.5, even 
with 8-meter class telescopes. 

When compared with the observational fit of 



IShaplev et alJ i200ll) . it is clear that the 'Q '-series are 
deficient in the brightest galaxies at the high luminosity- 
end of the luminosity function. This can be understood 
as a result of the small box size (ibox = 10 ft -1 Mpc) of 
these runs, which do not have large enough volume to 
allow a faithful sampling of rare, bright objects. As the 
box size becomes larger from the 'Q'-series to the 'D'- 
series (Lbox = 33.75 h~ x Mpc), and then to the 'G'-series 
(ibox = 100 h" 1 Mpc), this situation improves however. 
More and more of the luminous objects can then be found, 
and the agreement with the observation becomes better at 
the bright-end of the luminosity function. 

On the other extreme of the luminosity distribution, we 
see that increasing the resolution from Q3 to Q4, and then 
to Q5, allows inclusion of ever fainter objects, as expected. 
Therefore the luminosity function becomes wider towards 
the fainter end. Note, however, that the bright-end hardly 
changes, suggesting good convergence in the simulation re- 
sults for the massive galaxies. 

The run '03' (no wind run) slightly overpredicts the 
number of galaxies compared to observations, arguing for 
the need of stronger feedback. To investigate this further, we 
explicitly compare in Figure |S| the rest-frame V-band lumi- 
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Figure 5. Rest-frame V-band luminosity functions. The three different lines are for three different values of extinction: E(B — V) = 0.0 
(blue long-dashed), 0.15 (green solid), an d 0.3 (red dot-dashed ). The black short-dashed curve is the observational estimate of the 
rest-frame V-band luminosity function by IShaplev et all 1200 f}| . and the observational magnitude limit of My = —20.5 is shown by 
the arrows. Note that the Schechter function is plotted down to the fainter magnitudes only to guide the eye for comparison with the 
simulation results. 



nosity functions as a function of wind strength. We show the 
results for the simulations 03 (no wind), P3 (weak wind), 
and Q3 (strong wind), using an extinction of E(B— V) = 0.3. 
For reference, we i nclude the observational estimate by 
Shap levet al.l (l200ll) as a black short-dashed curve. Clearly, 
stronger galactic winds reduce the number of bright galaxies, 
as a result of the suppression of star formation they incur. 
The same effect has been found in the Eulerian hydro dy- 
namic simulations bv lNaeamine. Cen. fe Ostrikerl (120041 in 
preparation). 

6.2 Luminosity functions as observed at z=0 

In Figure|7| we show the .R-band luminosity function of sim- 
ulated z — 3 galaxies as seen at the present epoch. As in Fig- 
urc[5] the three different lines are for three different values of 
extinction: E(B — V) = 0.0 (blue long-dashed), 0.15 (green 
solid), and 0.3 (red dot-dashed). The open square symbols 
give the direct observational result for the luminosity func- 
tion of LBG galaxies in the observed frame, while the solid 
circles are the data corrected for dust extinction, as given in 
Fig. 13 of lAdelbereer fc Steidell feOOOl) . For the uncorrected 



data points (open squares), a Schechter function fit with 
parameters $* = 4.4 x 10~ 3 /i 3 Mpc~ 3 mag _1 , M* = 24.54, 
and a = —1.57 is also drawn as a short-dashed line. As 
discussed in Section |SJ we applied the lCalzetti et all ll200Cfl 
extinction law locally at z — 3 to the simulated galaxies, and 
then redshifted their spectra, while simultaneously account- 
in g for the absorp tion by the IGM following the prescription 
of iMadaul il995l) . The combined effects of redshifting and 
absorption by the IGM change the shape of the luminosity 
function in a complex manner compared to the rest-frame. 

For the Q-series, the agree ment between the simu l ations 
and the observational data of lAdelbereer fc Steidell <|2000T) 
is not very good. On the bright end, these runs simply lack 
galaxies due to their small box size, explaining the discrep- 
ancy. But the well resolved faint-end appears to show some 
excess over the observational result. The G-series agrees bet- 
ter with the observations however. In particular, the G5 
and G6 runs with E(B — V) = 0.15 agree well with the 
uncorrected observed luminosity function, and the case for 
E{B — V) = 0.0 is reasonably clo se to the dust corrected 
one oflAdelbereer fc Steidell feOOOD . 

lAdelbereer fc Steidell i200Q|) corrected their observed 
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Figure 6. Comparison of the rest-frame V-band luminosity func- 
tions of the 03 (no wind), P3 (weak wind), and Q3 (strong wind) 
simulations for E(B — V) = 0,3. The black short-da shed curve 
is the observational estimate by IShaplev et alJ 1200 ll) , and the 
observational magnitude limit of My = —20.5 is shown by the 
arrows. Note that the Schechter function is plotted down to the 
fainter magnitudes only to guide the eye for comparison with the 
simulation results. 

luminosity function for dust by adopting the relation 
A laoo = 4.43 + 1.99/3 ijMeurer. Heckman. fc OaIzettilll99flh 
for the extinction in magnitudes at 1600 A, where f3 is the 
UV slope of the spectrum which has a distribution in the 
range of [—3, 0] with a peak around —1.5, corresponding 
to_j4i6oo ~ 1-4 mag (see Fig. 12 of lAdelberger fc Steidell 
2000). This peak value is roughly consistent with the value 
of E(B — V) = 0.15 we adopted here with the lCalzetti et all 
( 2000 ) extinction law k(X), because E(B-V)k(W00A) ~ 1.5 
mag. Therefore the rough agreement between the simulated 
luminosity function with E(B — V) = 0.0 for G5 and G6-run 
and th e dust corrected data points of lAdelberger fe Steidell 
(2000) is encouraging, although not perfect. 

Also, note that the peaks of the luminosity functions of 
the G5 run with E(B — V) = 0.0 are on the brighter side 
of R = 25.5 owing to the lower resolution. This results in 
the increase of the number density of the LBGs that satisfy 
the colour-selection criteria when the resolution is increased 
from the G5 to G6 run as described in Section 2] 

Overall, the comparison of the simulated luminosity 
functions with the observational results presented in Sec- 
tions 16.11 and 16.11 stresses the importance of having a large 
simulation box size (> lOO/i" 1 Mpc) in order to obtain rea- 
sonable agreement at the bright end. In our simulations, only 
the G-series contain a sufficient sample of LBGs brighter 
than R = 25.5 with E(B - V) = 0.15. 



7 STAR FORMATION HISTORY OF LBGS 

In Figure|HJ we show four examples of typical star formation 
(SF) histories of galaxies in the 'G6' run with a bin-size 



of 10 Myrs, as derived from the age distribution of stars 
found in each galaxy including all the progenitors that merge 
prior to z = 3. We here chose the G6 run because it gives 
a reasonable agreement with observations both for the R- 
band and rest-frame V-band luminosity functions, and it has 
higher resolution than the G5 run. On the right hand side 
of each panel, we indicate for each galaxy an ID, its stellar 
mass in units of /i -1 M0, its apparent R magnitude (for 
E(B — V) = 0.15), and its rest-frame V-band magnitude. 

A notable feature present in all the panels is the numer- 
ous spikes of starbursts lying on top of a relatively continu- 
ous component. The starbursts last 10 — 20 Myrs, but the un- 
derlying smooth component shows only moderate variations 
in its rate when averaged over bins of 100 Myrs beginning 
from high redshift (z ~ 10) to z = 3 in these galaxies. For 
example, the galaxy shown in the top panel, which is one of 
the most massive galaxies in the G6 run with a stellar mass 
of M s tar = 1.4 x 10 11 h, _1 M Q and rest-frame V-band mag- 
nitude of —24.4, has continuously formed stars from z = 10 
to z = 3 at a typical rate of ~ 40 Mq yr _1 . 

In the other panels, we show galaxies of progressively 
smaller total mass, with correspondingly lower levels of star 
formation. Note that the galaxy on the bottom has a stellar 
mass of only M* = 3.2 x 10 10 /i~ 1 M Q , yet it still satisfies 
the colour-colour selection criteria, although its luminosity 
and star formation rate are substantially smaller compared 
with the galaxy shown in the top panel. The recent star 
formation near z — 3 has helped to bring the brightness of 
this relatively small mass galaxy above the limit of R = 25.5. 

Interestingly, the underlying continuous component 
of the star formation histories m easured here are much 
smoother than the ones found by iNagaminel ll2002ll when 
analysing a Eulerian hydrodynamic simulation. Their star- 
bursts were more extended on the order of 100 Myrs, and 
more sporadic. This difference is likely related to the differ- 
ent models used by the codes for the treatment of the physics 
of star formation. In the SPH methodology we investigate 
here, a sub-resolution multiphase model for the ISM was 
used which has a self-regulating property; i.e. gas that cools 
smoothly onto the ISM is also consumed in a smooth fash- 
ion by star formation, and only gas-rich major mergers can 
trigger starbursts. 



8 DISCUSSION 

We have used state-of-the-art hydrodynamic simulations 
of structure formation to study the properties of LBGs 
in a ACDM universe. Our simulations use a new 'entropy 
conserving' SPH formulation that minimises systematic in- 
accuracies in simulations with cooling, as well as an im- 
proved model for the treatment of the multiphase structure 
of the ISM in the context o f star formation and feedback 
jSpringel fc Hernaui st 2003a). For the first time, our study 
uses a large series of simulations to investigate the properties 
of LBGs, allowing an exploration of an unprecedented large 
dynamic range in both mass and spatial scales, while simul- 
taneously providing reliable estimates of systematic effects 
due to numerical resolution. 

By comparing our results for LBG properties in simu- 
lations with different box size, we showed that a sufficiently 
large volume (at least Lbox > 50 ft -1 Mpc) is crucial to faith- 
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Figure 7. R-band luminosity functions as observed at z = 0. The three different lines are for three different values of extinction: 
E(B — V) = 0.0 (blue long-dashed), 0.15 (green solid), and 0.3 (red dot-dashed). The ope n square symbols give the ob served luminosity 
function, while the solid circles are the dust corrected data points, as shown in Fig. 13 of Adelbcrecr & Stcidcl ( 20 



fully sample the LBG population. In particular, the bright- 
end of the LBG luminosity function is invariably incomplete 
in simulations with small box size. Of equal importance is 
a proper treatment of effective feedback processes. By con- 
sidering a series of runs with different strengths of galactic 
winds, we showed that without an effective feedback pro- 
cess, the number of LBGs is overpredicted. However, for our 
default model of strong winds, a reasonable space density of 
bright LBGs is obtained. In particular, we found that our 
G-series with a box size of Lbox = 100 h~ x Mpc has a quite 
plausible LBG population, with luminosity functions in both 
rest-frame V^-band and observed-frame it-band that match 
the observations reasonably well, at least at the bright-end. 

Perhaps the most important conclusion of this pa- 
per is that the observed properties of LBGs, including 
their number density, colours, and luminosity functions, 
can be well explained if the LBGs are simply associ- 
ated with the most massive galaxies at z = 3, with me- 
dian stellar mass of M, ~ 10 10 /i _1 M . This conclusion 
is consistent with earlier numerical studies b ased on hy- 
drodynamic simulations of the ACDM m odel (lPaveet_alJ 
| 1999(: iKatz. Hernquist fc Weinberg! Il999t iNagamind 120021: 
IWeinberg. Hernquist fc Katzl 120021) as well as some semi- 



analytic models teaueh et alJll99ctlKauffmann et al.ll999l) . 

and does not provide direct support for alternative models 
which suggest that LBGs are star-bursting low-mass sys- 
tems that later evolve into low-mass spheroids at z = 0. 
This point is corroborated by the high rates of star forma- 
tion with > 10 Mq yr _1 seen over extended periods of time 
of order 1 Gyr in the simulated galaxies, leading to the build- 
up of typical stellar masses of 10 10 h^MQ at z = 3. These 
comparatively steady star formation histories are also con - 
siste nt with observ a tional studies bv lPapovich et all (120011) 
and IShaplev et alJ (120011) . Note however that these same 
observations also find some evidence for starburst activity 
(> lOOMoyr -1 ) within the last 500 Myr before z = 3, last- 
ing for periods of ~ 100 Myr. Such long violent bursts are 
not seen in our SPH simul ations, but they w ere present in 
the simulations analysed bv lNagam inc (2002). Instead, mul- 
tiple shorter bursts with time-scale of ~ 10 — 20 Myr are 
seen in our SPH simulations. 

A recent study by Wea therlev fc Warrenl <l2003l) 
however cla i ms t hat the kinematic measurements of 
iPettini et all fcOOlj) and lErb et ail ll2003ft favour the picture 
of LBGs being low-mass starbursting systems. But the mea- 
sured 'rotation curves' are strongly dependent on seeing, and 
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Figure 8. Star formation histories of selected galaxies in the 
'G6'-run with the bin-size of 10 Myrs. Galaxy ID, stellar mass in 
units of /i _1 Mq, rest-frame V-band magnitude, and apparent R 
magnitude for E(B — V) = 0.15 are indicated in each panel for 
the corresponding galaxy. 



it is very likely that these measurements underestimate the 
true rotation curves (M. Pettini, private communication), 
casting some doubt on these claims. To settle this issue, we 
have to wait for the detection of a definite flattening of the 
rotation curves, which might become possible in the near fu- 
ture with near-IR spectrographs coupled to adaptive optics 
systems. Another difficulty with the low-mass starbursting- 
galaxy picture for LBGs is that such systems would be ex- 
pected to have relatively low chemical abundances, whereas 
at least the brighter LBGs have near-solar metallicity. 

The luminosity functions we measure for our simulated 
LBG galaxies are characterised by a steep faint end, with a 
slope close to a ~ — 2. This is the slope of the dark mat- 
ter halo mass function, so in our simulations, the galaxy 
luminosity function seems to largely follow the halo mass 
function at the low-mass end. Observationally, we know 
that at 2 = the faint-end slope of galaxies is certainly 
much shallower than this, with a slope of order a ~ —1.2 
(Blan ton et alJl2001l) . At high-redshift, the faint-end of the 
galaxy luminosity function is not nearly as we ll constrained, 
howev er, and may be much steeper. In fact, IShaplev et alJ 
tefJOll) report a slope of a = —1.85 for LBG galaxies at 
z — 3 which is still consistent with our simulations. How- 
ever, as the authors explain, this should be taken with some 
caution because their observations only reach an absolute 
1^-band magnitude of —20.5, which is not deep enough to 
constrain the faint-end slope reliably. If confirmed, a steep 
luminosity function at z ~ 3 would in any case have to 
evolve into a much flatter one at low redshift. Whether such 
a scenario is viable or not will be tested by current and fu- 
ture very deep surveys , such as DEEP 2 dPavis et alJl2002t 
iMadeewick et al.ll2003l : ICoil et aljEooll) . 



Barring confirmation of a rapid evolution of the faint- 
end slope, our results therefore hint that the simulations 
still overproduce the number of low-mass galaxies, despite 
the inclusion of strong feedback processes that are capable of 
accounting for the observed brightness of L*-galaxies. The 
moderating effects of feedback on star formation do depend 
on galaxy mass in the physical model followed by the simula- 
tions in the sense that, for a given amount of star formation, 
the feedback by winds does comparatively more damage for 
small mass galaxies. This is because towards smaller galaxy 
mass scales, the winds find it easier to escape from their 
confining galactic potential wells. The winds also entrain 
more gas in the outflow for lower circular velocities, so that 
the net baryonic loss becomes larger towards smaller mass 
scales, making the total mass-to-light ratio of smaller halos 
in principle bigger. But, this variation of the feedback effi- 
ciency with galaxy mass may not be strong enough in the 
present set of simulations to close the significant gap be- 
tween the faint-end slopes of halo mass function and galaxy 

luminosity function. 

In terestingly, the s t udies by IChiu 1 Gnedi n fc Ostrikerl 

j200ll) : lNagamine et all feOOll) . and lNaeamind d2002l) based 
on Eulerian hydrodyn amic simulations did fin d a flatter 
faint-end slope (see also lHarford fc Gnedinl20 03^ . Presently 
it is unclear whether this was due to resolution limitations 
in the low-mass end of the halo mass function, or due to 
genuine physical effects of the feedback model implemented 
in these Eulerian simulations. More work in the future will 
be needed to settle this very interesting question, which is 
of tremendous importance for the theoretical framework of 
hierarchical galaxy formation in CDM universes. 
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